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ABSTRACT 

We present a new method for constructing equilibrium phase models for stellar sys- 
tems, which we call the iterative method. It relies on constrained, or guided evolu- 
tion, so that the equilibrium solution has a number of desired parameters and/or 
constraints. This method is very powerful, to a large extent due to its simplicity. It 
can be used for mass distributions with an arbitrary geometry and a large variety of 
kinematical constraints. We present several examples illustrating it. Applications of 
this method include the creation of initial conditions for TV-body simulations and the 
modelling of galaxies from their photometric and kinematic observations. 

Key words: galaxies: kinematics and dynamics - methods: N-body simulations 



1 INTRODUCTION 

In astronomy there are at least two problems where equilib- 
rium phase models of stellar systems need to be constructed. 
The first one is the construction of phase models for real 
galaxies from observational data, i.e. the modelling of ob- 
servational data. The second problem is the construction of 
initial conditions for A-body simulations of stellar systems. 
It is obvious that these two problems are tightly connected, 
and yet they have, so far, been solved by different m eth- 
ods. The Schwarzschild method (|Schwarzschildlll979t ) and 
its modific ations is often used f o r modelling of observational 
data (e . g lHafner et all l2000l: Ivan den Bosch et all l200rj; 
Thomas! 120071 " Ivan den Bosch et al.l l200Sl : Ide Lorenzi et al.l 
20081 '). but has almost never been used so far to produce 



initial conditions for simulations. For .A-body initial condi- 
tions, a wide variety of methods has been used, based on the 
Jeans theorem (e.g. IZangjll976l ; lAthanassoula fc Sellwoodl 
19861: iKuiiken fc Dubinskill 19951 : Iwidrow fc Dubinskill2005l ; 
' 120071 ) . 



McMillan & D ehnenl |2007|), or on Jeans' equations (e.g. 



HernquistJll993l ). In the case of multi-component systems, 
e.g. disc galaxies with a bulge and a halo, the components 
are built separately and then either simply superposed, or 
the p otential of the one is adiabatically grown in the other 
(e.g iBarnesll 19881 : iMcMillan fc Dehnenll2007l : lAthanassoula! 
120071 ) before superposition. 

For real galaxies the phase space density is generally 
unknown, but we do have some information about it. For 
example, we know more or less accurately a distribution of 
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mass for the visible components (notwithstanding uncertain- 
ties due to the mass to light ratio) and we often have some 
constraints on the velocity distribution. It is also reason- 
able to assume that the galaxy is in an equilibrium state. So 
in general, the problem of constructing a model in phase 
space is equivalent to constructing an equilibrium phase 
model with a given mass distribution and, in many cases, 
given kinematic constraints. In the case of modelling obser- 
vational data (first of the two above mentioned problems) 
the kinematic parameters are the observed velocities inte- 
grated along the line of sight. In the case of initial conditions 
for iV-body simulations, a wide variety of kinematic param- 
eters is possible, depending on the problem the simulation 
addresses. 

We have developed a new method for constructing equi- 
librium phase models with a given mass distribution and 
with given kinematic parameters, which we call the iterative 
method. It can be applied to systems with arbitrary geom- 
etry, so that the requested mass distribution can be arbi- 
trary. The idea and a first im ple mentation was prese n ted in 
iRodionov fc Sotnikoval (|2006u . In lRodionov fc Orlovl (I2008D 
we improved it, and applied it to construct an A-body model 
of the stellar disk of our Galaxy for two realistic mass mod- 
els of the Milky Way. Here we present a final version of this 
method, fully allowing kinematical constraints. In the previ- 
ous articles we had concentrated on constructing equilibrium 
phase models with a given mass distribution, so that kine- 
matic parameters were either not considered or only in terms 
of au xiliary parameters, such as the total angular momen - 
tum (|Rodionov fc Sotnikoval2006l : IRodionov fc Orlovll2008l ) . 
This, however, limited the applicability of our method, both 
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for initial conditions and for modelling real galaxies. Indeed, 
initial kinematics play a crucial role in determining the evo- 
lution of TV-body systems, while observational constraints 
more often than not include kinematics. In this paper we 
give equal attention to the mass distribution and kinemat- 
ical constraints, so that the iterative method can now be 
used for a number of interesting applications. In principle, 
in our method, both the kinematic constraints and the mass 
distribution can be arbitrary. But the part of our algorithm 
that concerns the kinematic constraints is not universal, con- 
trary to the part that handles the mass distribution, but is 
tailored to the specific constraint. Here we consider several 
types of constraints. Once these are understood, it is rather 
easy to extend the algorithm for every new type of kinematic 
parameter (see below). 

The power of the iterative method stems from its sim- 
plicity. The iterative method is based on a simple and, in a 
way, obvious idea, which is implemented in a simple algo- 
rithm. The purpose of this article is to fully describe this 
method. We first introduce the basic concept in Section [2] 
where we also explain the different modules of the algorithm 
and the way they should be applied. In section [3] we illus- 
trate the use of the method with three examples, namely a 
triaxial system, a multi-component model of a disk galaxy 
(including live disk, bulge and halo components) and a disk 
constructed with given line-of-sight kinematic. We briefly 
conclude in section [4] 



2 THE ITERATIVE METHOD 

2.1 General Idea of the Iterative Method 

The aim of the iterative method is to construct equilibrium 
TV-body models with a given mass distribution and with 
given kinematic properties, parameters, or constraints. This 
method relies on the fact that any non-equilibrium system 
will tend, more or less fast, to a stable equilibrium. We thus 
start by constructing any arbitrary, non-equilibrium TV-body 
system, and let it evolve. Such an evolution changes both its 
mass distribution and its kinematics, so that the final sys- 
tem does not have the desired properties. To achieve the 
latter, we developed a new method which we call the iter- 
ative method and which relies on a constrained, or guided, 
evolution. We will describe it fully in this section. This idea 
is of course applicable for any arbitrary dynamical system 
and is even widely used in every day life. 

To give an example, let us consider a donkey walking by 
itself in a field. After some time the donkey can be anywhere 
in the field. Now consider another donkey which we attach 
to a tree by a rope. This donkey also walks in the field, but it 
will have to stay inside a circle of radius equal to the length 
of the rope. This is an example of a constrained evolution, 
which will necessarily lead to a final state within a circle 
around the tree. The crucial point now is how to achieve 
this constraint, i.e. what will the equivalent of the rope be 
in the case of galaxy models. 

The general scheme of our method is presented in Fig. [T] 
and can be applied to any arbitrary dynamical system. Let 
our task be to find an equilibrium state of some dynamical 
system obeying given constraints or having specific values of 
some given parameters. We start from any arbitrary state of 
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Figure 1. General scheme of the iterative method in the case of 
an arbitrary dynamical system. 
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Figure 2. The scheme of the iterative method for the case of an 
TV-body system with a given mass distribution and given kine- 
matical parameters. 

our dynamical system and allow the system to evolve dur- 
ing a short time interval. We then need to make sure that 
the given parameters have the required values. In order to 
achieve this we need to modify the system so that the given 
parameters have the necessary values, while making sure 
that the other quantities or parameters are kept unchanged, 
so as to retain memory of the evolution. As shown in Fig. [T] 
we iterate these two steps, alternating short evolutions and 
modifications of the system to fix the set parameters. We 
thus constrain the evolution in order for it to reach an equi- 
librium state with the desired set of constraints. We stop 
when we consider we are sufficiently near the desired equi- 
librium state of the system. 

Let us now consider a case, in which we wish to con- 
struct an equilibrium TV-body system with a given mass dis- 
tribution and with or without given kinematic constraints. 
The scheme is outlined schematically in Fig. [2] We initially 
create an TV-body system with a given mass distribution 
but with arbitrary particle velocities (for example veloci- 
ties equal to zero). We then start the iterative procedure, 
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by letting the system go through a sequence of evolution- 
ary steps of short duration. At the end of each one of these 
steps, and before the new step is started, we need to set 
the appropriate parameters. Let us first consider the case 
where we wish to have a specific mass distribution, but have 
no kinematic constraints. To achieve this, we construct a 
new TV-body system, with the desired mass distribution but 
with velocities chosen according to the velocity distribution 
obtained from the evolution. In other words, we "transfer" 
the velocity distribution from the system obtained from the 
evolution to a new system, which will have the desired mass 
distribution and an evolved velocity distribution. The algo- 
rithm performing this "transfer" is the core of the iterative 
method, and will be discussed in more detail in section [2~2l If 
we have kinematic constraints as well, we need to modify the 
velocities of the particles in such a way that the constraints 
are fulfilled and, at the same time, as little as possible so 
that some memory of the evolution is kept. How this is done 
in practise depends on the imposed constraints and will be 
described, for a number of cases, in section [2731 Procedures 
for further types of constraints can be easily found follow- 
ing similar techniques. In all cases we have a new system 
which has the desired mass distribution, obeys the neces- 
sary velocity constraints, while being nearer to equilibrium, 
since it retains partial memory of the evolution. We repeat 
this iterational procedure a number of times, alternating one 
evolution phase and one phase where the necessary parame- 
ters are set, until we come as near to the desired equilibrium 
state as desired. 

2.2 Transfer of velocity distribution 

The transfer problem can be formulated as follows. Any evo- 
lution step ends with a model, which we will refer to as the 
"old" model. This step is followed by a constraining, or fixing 
step, during which we create a "new" model with the desired 
mass distribution. We now need to transfer the velocity in- 
formation from the "old" to the "ne w" model. There is more 
than o ne way to achieve this transfer. iRodionov fc Sotnikoval 
(2006) used an algorithm based on moments of the velocity 
distribution, which, however, proved to be rather compli- 
cated and cumbersome. Here we suggest a much simpler and 
more reliable algorithm, whic h is in fact an impr o veme nt of 
an algorithm initially used in lRodionov fc Orlovl (|2008T ). 

The basic idea of our velocity transfer algorithm is as 
follows. We assign to the new-model particles the veloci- 
ties of those particles from the old model that are nearest 
to the ones in the new model. The simplest (although, as 
we show below, not necessarily optimum) implementation 
of this idea is evident. One can prescribe to each particle in 
the new model the velocity of the nearest particle from the 
old model. Let us formulate this proposition more strictly. 
For each i-th particle of the new model, one finds the j- 
th particle in the old model with the minimum value of 
\rf ew - rf d \. Here, r" e ™ is the radius vector of the i-th par- 
ticle in the new model, and r° l is the radius vector of the 
j-th particle in the old model. Hereafter we will always imply 
this definition when we talk of the nearest or closest parti- 
cle. One then takes as the velocity of the i-th particle in the 
new model the velocity of the j-th particle in the old model. 
This simple algorithm has, however, one essential drawback. 
If the numbers of particles in the old and new models are 



the same then only about one-half of the particles in the old 
model participate in the velocity transfer. The reason is that 
many old-model particles transfer their velocities to several 
particles in the new model. As a result of this, almost one- 
half of the particles in the old model do not transfer their 
velocities at all. This means that a significant amount of 
information on the velocity distribution will be lost in the 
transfer process. The noise will therefore grow, as we veri- 
fied in numerical experiments. Thus, this transfer algorithm 
is not optimum. 

This shortcoming can, nevertheless, be overcome by 
modifying this transfer scheme. For this, we introduce an 
input parameter, which we call the "number of neighbours" 
n n b- We also introduce, for each particle in the old model, 
the parameter n use d, which denotes the number of times this 
particle has been used for velocity copying. At the beginning 
of the transfer procedure we set n uae d = for each parti- 
cle in the old model, since its velocity information has not 
been yet transferred to any of the new model particles. For 
any given particle in the new model we find the nearest n„b 
neighbours in the old model (the closeness being understood 
as defined above) and from these we single out the subgroup 
of particles that have a minimum n U sed- From this subgroup 
we find the particle that is the closest one to the position 
of the new-model particle, add one to its n use d value and 
prescribe its velocity to the new-model particle we are ex- 
amining. 

We note that if n n b = 1 this algorithm is the same as 
the previous one and about half of the particles will not take 
part in the velocity distribution transfer. If, however, we take 
n n b = 10, only a small fraction (a few per cent) of old-model 
particles will not take part in the transfer process. We adopt 
this improved transfer method since we showed that it gives 
good results in the iterative procedure. 

If the desired model has some symmetry, it can be useful 
to take it into account in the algorithm of velocity transfer by 
redefining the distance between two particles. For example, 
if we wish to build an axisymmetric system, we search for the 
nearest particles in the two-dimensional space R — z (where 
R the cylindrical radius) instead of the three-dimensional 
space x — y — z. We then transfer the velocity of this nearest 
old-model particle (in cylindrical coordinates) to the new- 
model particle. It is important to adopt this new definition of 
the distance in order to fix not only the mass distribution, 
but also fix the velocity distribution and to make it fully 
axisymmetric. 

We have thus introduced three variants of the velocity 
transfer algorithm. We will refer to them as "transvel_3d" , 
"transveLcyl" and "trasveLsph" . 

(i) "transvel_3d" : This is the basic algorithm, for the case 
when the desired system has no assumed symmetry. By us- 
ing this algorithm in the iterative method we only fix the 
mass distribution and leave the velocity distribution un- 
changed. In the current work we use this algorithm when 
constructing triaxial models. 

(ii) "trasnveLcyl" : This is a modification of the basic al- 
gorithm for axisymmetric systems and was described just 
above. We use this algorithm when both the desired mass 
and velocity distribution are axisymmetric. In the current 
work we use this algorithm for constructing all models ex- 
cept for the triaxial ones. 
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(iii) "transveLsph" : This is a modification of the basic 
algorithm for spherical systems. In this version of the algo- 
rithm we search for the nearest particles in one dimensional 
"r" space, where r is the spherical radius. By using this al- 
gorithm in the iterative method we fix both the mass and 
the velocity distribution to m ake the system fully sphericall y 
symmetric. This was used in ISotnikova fc Rodionovl |2008). 



2.3 Fixing the Kinematic Parameters 

Here we describe algorithms for fixing different kinematic 
parameters. The general algorithm is as follows. We slightly 
change the particle velocities to fix given kinematic param- 
eters, but we keep as many other parameters as possible 
unchanged. Here we describe in detail only a number of algo- 
rithms, which we use in this paper. But it is easy to develop 
algorithms for any other kinematic parameter. It is only nec- 
essary to follow the general principle: "keep unchanged the 
parameters that do not need to be fixed" . 



2.3.1 Fixing the radial velocity dispersion profile or(R) 

We use this algorithm in order to fix in stellar disks the 
radial velocity dispersion profile to a given function or(R) 
(for an application, see section ET2")) . 

Let o~fi(R) be the given radial velocity dispersion profile 
which we want to fix, where R is the cylindrical radius. Af- 
ter each evolutionary step (see Sect l2.ip we need to change 
slightly the radial velocities of particles in order to fix this 
profile. The model is divided into ridiv concentric cylindrical 
annuli, each containing the same number of particles. For 
each annulus j, we calculate the target value of the radial 
velocity dispersion 



o R = o R (Rj 



(1) 



where Rj is the mean value of the R coordinate of all par- 
ticles in part j. The new radial velocity of the i-th particles 
in the j region is then prescribed as follows. 



VRi 



I j j j/ 



(2) 



where v' Ri is the current value of the i-th particle radial 
velocity, vm is the corrected i-th particle radial velocity and 
o R is the current value of radial velocity dispersion in part 
j. We note that in this scheme we have assumed that the 
mean radial velocity is equal to zero. 



2.3.2 Fixing the radial anisotropy profile 

This algorithm is very similar to previous one and is very 
useful for building spherical models with a given profile 
of velocity anisotropy. We will use it for building the 
3.21 and it has also been used in 
20081). 



halo model in section 
ISotnikova fc Rodionovl 

Let ag, "V and ay be the velocity dispersions in the 9, 
if and r directions of spherical coordinate system and let us 
aim e.g. for a model with a given profile, P(r), of the ae/cy 
ratio. The model is divided in concentric spherical shells, 
each containing the same number of particles. For each shell 
j, we calculate the target value of ae/ay 



where rj is the mean value of the r coordinate of all particles 
in shell j. We will attempt to obtain this ratio by changing 
appropriately the 9 component of particle velocities (alter- 
natively, we could have changed the r component). The new 
6 velocity component of the ith particle in the jih region 
will then be prescribed by 



(4) 



where v' gi and voi are, respectively, the current and the cor- 
rected values of the i-th particle 9 velocity component and 
a' r and a'g are the current values of the radial and 9 velocity 
dispersion, respectively, in part j. 



2.3.3 Fixing the line-of-sight mean velocity or the 

line-of-sight velocity dispersion in the case of an 
edge-on disk 

In this section, we describe two algorithms, one for fixing 
the line-of-sight mean velocity and the other for fixing the 
line-of-sight dispersion of an axisymmetric disk. To set the 
notation, let us assume that the stellar disk rotates about 
the 2-axis, the disk plane lies in the {x, y) plane and the line 
of sight is along the y axis (edge-on disk). We invert the 
sign of v y for each particle with x < in order to make the 
half disk with x < kinematically identical to the half with 
x > and flip the x < particles on the x > part. We 
then divide the disk in slits parallel to the (y, z) plane and 
at different distances from the centre, i.e. at different values 
of x, in such a way that all slits have the same number of 
particles. 

Let us denote by v\ os (x) the desired profile of the line-of- 
sight mean velocity, i.e. the mean value of v y after integra- 
tion along the line of sight. For each slit j we calculate the 



target value of the line-of-sight mean velocity vf c 



Vl QS (Xj) 



(where Xj is the mean value of \x\ for particles in part j) 
and the current value of the line-of-sight mean velocity v^ s 
(as the mean value of v y for all particles in slice j). The new 
y velocity component of particle i in region j should then be 



3 . 
los 



(5) 



where v' yi is the current value of the i-th particle y veloc- 
ity and v y i is the corrected i-th particle y velocity. Particles 
which were flipped to x > part have to be flipped back and 
the sign of their y velocity component reversed. The particles 
are then azimuthally mixed to make the velocity distribu- 
tion axisymmetric and the step is concluded. Of course, in 
this way we have tampered with vi os , but this unavoidable. 
Nevertheless, after a number of iterations, both the axisym- 
metry and the desired vi oa will be achieved. 

The algorithm for fixing <ri os (a;) is very similar, except 
that we have to calculate in each slit the current value of the 
line-of-sight velocity dispersion cr? os as the dispersion of v y 
for all particles in slit j. Let a\ os (x) be the desired profile of 
the line-of-sight velocity dispersion. In order to achieve this, 
the new y velocities should be modified as follows 



{v'yi 



j, + U los 



(6) 



ft = P(rj 



(3) 



where of = a\ os (xj) is the target value of the line-of-sight 
velocity dispersion. We note that in this algorithm we have 
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to take into account that the value of need not neces- 

los 

sarily be equal to zero. As previously, we still have to flip 
back particles which were flipped to the x > part, invert 
the sign of their y velocity component with and mix the 
particles azimuthally. 

2.3.4 Fixing velocity isotropy conditions 

An isotropic velocity distribution depends only on the veloc- 
ity module and not on the direction of the velocity. Our algo- 
rithm for fixing it is very simple. For each particle, we keep 
the velocity module unchanged and randomise its direction, 
thus ensuring that the velocity distribution is isotropic. For 
spherical isotropic models, the distribution function (DF) is 
known, at least formally, or numerically. Thus the construc- 
tion of such models can be considered as a test of the itera- 
tive method, and we have verified in a number of cases that 
the models constructed by the iterative method are identical 
to the models constructed by using known equilibrium DF. 
The construction of spherical isotrop ic models with the iter- 
ative method was first described in iRodionov fc Sotnikoval 
(2006). That work, however, used an old and rather com- 
plicated algorithm for transferring the velocity distribution. 
Here we use a different, superior algorithm, based on the de- 
scription in Sect. [2721 We again made sure that the models 
thus constructed are identical to the models obtained by us- 
ing a known DF. Furthermore, we also used this algorithm 
for constructing models which are not fully isotropic models, 
but rather not-very-far from isotropic (see section 13.11 and 
GL2J. 

2.4 How many parameters should we fix? 

The goal of the iterative method is to construct equilibrium 
iV-body models with given parameters (i.e. with a given 
mass distribution and with given kinematic constraints). 
There are in general three possible cases with respect to 
the number of constraints. 

In the first case the number of given parameters, or con- 
straints are such that only one equilibrium model can exist. 
In this case, we can expect that the iterative method will 
converge to this unique equilibrium model, independent of 
the initial state from which the iteration is started. For ex- 
ample, it is known that for a spherical model with a given 
mass distribution only one isotropic equilibrium DF exists. If 
we construct a spherical model by using the iterative method 
and we fix velocity isotropy as a kinematical constraint (sec- 
tion [2]3]4]), then the iteration always converges to the same 
model, independent of the initial state, as expected. 

In the second case, the number of give parameters is 
such that many equilibrium models can exist, i.e. this num- 
ber is insufficient for determining uniquely the equilibrium 
model. In this case we can expect that the result of the iter- 
ative method will depend on the choice of the initial model. 
The iteration will converge to the equilibrium model which 
is "nearest" in some sense to the initial model. Alternatively, 
the iteration will converge to some specific, in some sense, 
model. For example, when constructing a triaxial model in 
section 13.11 we fix only the mass distribution and do not 
set any kinematic constraints. Of course in this case the re- 
sult of the iterations will depend on the initial state (see 



section [3~T1 for details). Another, more involved, example is 
the construction of a disk model with given total angular 
momentum. In principle, man y such equilibrium models are 
possible, yet the iterations of IRodionov fc Orlovl |2008l ) al- 
ways converged to the same model. It is unclear why this is 
the case, but it could be du e to a specificity of the model 
f see IRodionov fc Orlovll2008l for details). 

The last possibility is that for the adopted parameters, 
no equilibrium model exists, i.e. we have fixed too many 
parameters. In this case the iteration will either not converge 
at all, or it will converge to a system with the parameters 
we have fixed, which is in non-equilibrium, but close in some 
sense to equilibrium. 

2.5 Technical comments 

In this section we elaborate a few important technical points, 
useful for anybody wishing to apply the iteration method. 

One of the free parameters of the iterative method is the 
duration ti of each iteration, i.e. the time interval over which 
the system is evolved during each iteration. How should the 
numerical value of ti be chosen? It is clear that this time 
should not be too short, so as to allow the system to evolve 
sufficiently during one iteration step. On the other hand, it 
should not be too long either, so as not to permit the evo- 
lution of various instabilities; otherwise, these instabilities 
may change the system substantially. For example, when 
constructing a disk system it is necessary to use iteration 
steps considerably shorter that the growth time of the bar 
instability, which of course varies strongly from one model 
to another. For this reason, there is no strict criterion and ti 
should be chosen empirically. Our experiments have shown 
that it is usually better to try relatively big ti values, thus 
ensuring a much faster convergence. Moreover, in some sit- 
uations the iterations for relatively small ti don't converge 
at all, while iterations for relatively big ti do. This was, for 
example, the case when we constructed a model with rela- 
tively cold stellar disk. So if iterations don't converge or they 
converge too slowly, it is often useful to consider bigger ti 
(within of course reasonable limits). Examples of appropri- 
ate ti values are given in all examples in Sect. [3] Moreover, 
our simulations have shown that, if we take ti within reason- 
able limits, the result of the iteration is the same (within the 
noise limits) and independent of ti, provided of course the 
chosen number of parameters and constraints allow a single 
solution. If the latter is not the case, and the result of the 
iteration depends on its starting model, then of course the 
result can depend also on ti. 

Another parameter of the iterative method is the pa- 
rameter n„b in the algorithm of velocity transfer (see sec- 
tion [23]) • Its value has been chosen in a more or less "ad 
hoc" manner and, by trial and error, we have found that 
a value of n n t = 10 is often satisfactory. Our test simula- 
tions have shown that the results of the iterative method for 
n n i, = 10 and n n j, = 100 are practically the same, at least 
for a total number of particles as those used in our trials, 
i.e. of the order of a few hundred thousands to a couple of 
millions. 

The most computer costly part of our method is the 
computing of the evolution of the system in each iteration, 
since the computing cost of all other parts of the method 
is very small. For this reason it is recommended to use a 
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fast TV-bod y code and we have adopted gyrfalcON l|Dehnenl 
l200d , 120021 ). Furthermore, our test simulations have shown 
that computation of the evolution can be carried out with 
relatively low accuracy. This is mainly due to the fact that 
we need to calculate the evolution only during short time 
intervals, so that errors do not have sufficient time to accu- 
mulate. Therefore, in the iterative method we use Gyrfal- 
cON with relatively big values of the tolerance parameter 
9t and of the time-step (see section |3J). Usually the total 
computing cost is considerably smaller, but still of the same 
order as that necessary to run a simulation with the con- 
structed model. This of course will depend on whether we 
can start the iterations form a model reasonable close to the 
final one, or whether lack of any a priori knowledge leads 
us to start, e.g. from zero initial velocities. A 'trick' which 
helps reducing the computing cost is to make a few itera- 
tions initially with a small number TV of particles and then 
gradually increase TV to the required number. In the proce- 
dure of velocity transfer described in section [2~2l the number 
of particles in the old and the new system can be different. 
So in the next iteration step we can get a system with a 
larger number of particles. 



3 EXAMPLES OF MODELS 

In this section we consider three examples of models con- 
structed by our method, namely a triaxial model of a 
spheroid, a multi-component model of a disk galaxy and 
a model with given line-of-sight kinematics. 



3.1 Triaxial model 

Our first example is that of a triaxial model. As mass dis- 
tribution we use a Plummer sphere flattened in two dimen- 
sions. 

-5/2 

, (7) 



, , 3M p i 2 / x 2 y 1 z 1 2 
P(x, y, z) = - — dpi -j + tj + -j + a p i 



where M p i and a p i are the total mass and the scale-length of 
the model and a, b, c are rescaling parameters. In the present 
specific example we discuss a model with the following pa- 
rameters: M p i = 1, a p i = 37r/16, a — 1, b — 0.8 and c = 0.7. 
Scaling our model to an elliptical galaxy with a — 3 kpc and 
M p i = 10 11 M a , gives a time unit t u » 17 Myr. 

Our target is to construct an equilibrium TV-body sys- 
tem with this given mass distribution. We didn't impose any 
well-defined restriction on the kinematics of the system, but 
aimed instead for a velocity distribution not very far from 
isotropic. Our initial model was cold with velocities equal 
to zero. Our target model was reached in 50 iterations. For 
the first 10 we fixed a condition of velocity isotropy (sec- 
tion [2]3]4}, while for the last 40 we didn't fix any kinematic 
parameters. If we performed all 50 iterations without fixing 
any kinematic parameters, we would also obtain an equi- 
librium model but with a higher degree of anisotropy. We 
chose TV = 500 000 and an iteration time fej = 10. The inte- 
gration step and the softening length were taken dt — 1/2 7 
and e = 0.01, respectively. The tolerance parameter for gyr- 
falcON was set to 9 t = 0.9 (see section 12 . 5|) and we used 
the "transvel_3d" modification of the algorithm of velocity 
transfer (see section [2~2)l . 
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Figure 3. Evolution of the ellipticity for the three projections 
of the triaxial model constructed with our iterative method. Note 
that they do not evolve, i.e. that the model is in equilibrium. 



In such a case it is crucial that the final model be suf- 
ficiently close to equilibrium so that the axial ratios do not 
tend to unity after some time evolution, as it happens for 
many other techniques used in calculating triaxial equilib- 
ria. To test this we evolved our model for 50 time units, i.e. 
roughly 50 crossing times for the scale-length of the system 
dpi. The integration step and softening length were taken 
dt — 1/2 8 and e = 0.00 5, respectively - in agreement with 
the recommendations of Rodionov fc Sotnikova! l|2005h (see 
also lAthanassoula et al.l2000r ) - and the tolerance parameter 
for gyrfalcON was set to 6 t = 0.6. Fig.[3]shows the evolution 
of the ellipticity of the model for three different projections. 
This was calculated as the ratio of the medians of the abso- 
lute values of corresponding particle coordinates. As can be 
seen, the shape of the model is practically unchanged dur- 
ing the evolution. We also made sure that the model also 
conserved all its other properties, thus demonstrating that 
it is indeed very close to equilibrium. 



3.2 Multi-Component model of a disk galaxy 

As a second example, we constructed a model of a disk 
galaxy consisting of three components: a stellar disk with 
a given profile of radial velocity dispersion, a non-spherical 
bulge and a halo with a given anisotropy profile. 

To start, we need to define the mass distribution in each 
of the components. The disk model is an exponential disk 
with density: 



Pd(R,z) = 



■ exp 



(-£) sech2 (£)' 



(8) 



where Md is the total disk mass, Rd is the disk scale length, 
Zd is its scale height and R is the c ylindrical radius. Th e 
halo model is a truncated NFW halo l|Navarro et al.lll996t ) 



Ph(r) = C h 



exp(-r 2 /r 2 h ) 
(r/r h )(l + r/r h y 



(9) 



where rn is the halo scale length, Ch is a parameter defining 
the mass of the halo and r t h is the truncation radius of 
the halo. For the bulge we used a truncated and flattened 
Hernquist sphere l|Hernquistl lT990') with density 
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Figure 4. Rotation curve for the disk galaxy model. The solid 
line is the total rotation curve. We also show the contributions 
from the disk, halo and bulge components. 




Figure 5. Radial profile of the Toomre parameter Q for the disk 
in the model described in Sect. 13.21 



Pb(R,z) 



where 



M' h r h exp{-d 2 /rl h 
2irq d(d + r b ) 3 



d = 



R2 + ^> 



(10) 



(11) 



r b is the bulge scale length, Ml is the total bulge mass before 
truncation, r t b is the truncation radius of the bulge and q is 
the flattening parameter. 

For the parameter values we chose for the disk: Md = 1, 
R d = 1, z = 0.2, for the halo: r h = 4, C h = 0.01, r th = 14 
and for the bulge: r b = 0.2, M' b = 0.2, r tb = 2, q = 0.7. 
For these parameters the total mass of the bulge and of the 
halo are Mb ~ 0.15 and Mu ~ 4.98, respectively. We use 
units such that the constant of gravity is G — 1. Scaling 
our model to a disk galaxy with Rd = 3.5 kpc and Md = 
5 ■ 10 10 Mq, gives a time unit t u w 13.8 Myr and a velocity 
unit v u ps 247.9 km/s. The rotation curve for our model is 
shown in Figure [3] which also displays the contribution of 
the disk, halo and bulge components separately. 

To make the exercise more realistic, we still need to 



choose kinematical constraints for each of the components, 
although, as we have already mentioned, these are not oblig- 
atory for our method. We created the disk with the following 
profile 



ctr(R) = 0.3 ■ exp (-R/3) + 0.2 ■ exp (-R/Q.5) 



(12) 



where gr is the radial velocity dispersion. From the mass 
model of the galaxy and from profile (| 12f) w e can calculat e 
the radial profile of the Toomre parameter Q (|Toomrelll964r ) , 
shown in figure [S] Our main target here is to demonstrate 
that our method can construct an equilibrium model of the 
disk with any realistic profile of <7r(R). The choice of profile 
in eq. (|12|1 is more or less arbitrary, but demonstrates that 
our method can work with more elaborate profiles than a 
single exponential function. 

When constructing the bulge, we did not impose any 
specific kinematic constraints. Instead, we aimed for a model 
not- very- far from isotropic, as in the case of the triaxial 
model of Sect. 13.11 

For the halo we adopted a velocity anisotropy profile, 
so as to test a different kind of constraint. More specifically 
we chose 
oeir) 



0.2 



.(r) 



+ 0.8, 



(13) 



r 

079 



^ +1 



where as and oy are the velocity dispersion in the 9 and 
the r direction in a spherical coordinate system. Note that 
we constructed the phase model of the halo in the presence 
of the non-spherical potential of the disk and bulge, i.e. we 
don't have spherical symmetry and in the halo equilibrium 
model a$ should not be equal to a v (the velocity disper- 
sion in <p direction) . So when we constructed the halo model 
we fixed only the fraction ag/a r , but did not fix the frac- 
tion o v /a r . This is different from t he case of the isolated 
spheri cal NFW halo, constructed in lSotnikova fc Rodionovl 
(2008). And again we want to underline that kinematical 
constraints are not obligatory. We could also construct the 
halo, or the bulge, without any kinematical constraints, or 
with another type of kinematical constraints. For example 
we can construct a rotating halo with given total angular 
momentum or a rotating halo with a given profile of the 
mean azimuthal velocity. 

Once the mass models and the required kinematical pa- 
rameters for each of three components are defined, we can 
apply the iterative method for constructing an equilibrium 
iV-body model for the whole system. In this specific exam- 
ple we chose N d = 200000, N b = 30324, N h = 995978 for 
the number of particles in disk, bulge and halo, respectively. 
With these numbers, the mass of particles in all components 
is the same. We constructed each of these components sepa- 
rately in the rigid potential of all other components. In order 
to take into account the external potential, we need to make 
only one small evident modification of the iterative method. 
Namely, when we need calculate the evolution of the system 
during the iteration time (see fig. [2} , we simply need to do it 
in the presence of the external potential. This can be done 
either by introducing an analytical external potential to the 
gyrfalcON program, or we can add it as a rigid iV-body 
system. In current work we use the latter. For example, in 
order to add a rigid halo we simply add in the system rigid 
particles according to the mass distribution of the halo. 
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Let us first describe the disk construction. Our initial 
model was a cold disk where all particles move on circu- 
lar orbits. Indeed, the circular velocity can be easily calcu- 
lated, since the mass distribution in the model is known. 
Had we, instead, started off with zero disk velocities, we 
would have again obtained an equilibrium model, but with 
counter rotating subsystems. This will happen because we 
fix only the profile of <jr(R) and do not fix any parameters 
defining the direction of rotation. It is therefore better to 
start off the iterations with a rotating disk, unless of course a 
disk with counter-rotating components is specifically sought. 
Note that the result of the iteration is independent of the 
initial iterative guess for the disk rotation. For example it 
will be the same if initially all the disk particles have tangen- 
tial velocities equal to half of the circular velocity, or twice 
that. 

We made 50 iterations, each with ti = 20. The inte- 
gration step and softening length were taken dt = 1/2 4 and 
e = 0.04, respectively and the tolerance parameter for gyr- 
falcON was set dt = 0.9. In order to fix the <jr(R) profile we 
used the algorithm described in section 12.3.11 The number 
of layers in this algorithm was n^iv = 200 (see section l2.3.ip . 
We used the "transveLcyl" modification of the algorithm of 
velocity transfer (see section 12. 2\ . This algorithm also was 
used for constructing the bulge and halo components. We 
call this disk model DISK.SVR. 

For constructing the bulge we also made 50 iterations. 
The other parameters for this construction were ti — 10, 
dt = 1/2 6 , e = 0.02 and 6 t = 0.9. Our initial model was 
a cold model with velocities equal to zero. During the first 
10 iterations we fixed a condition of velocity isotropy (sec- 
tion [2]3]4]), while we did not set any kinematical constraints 
during the last 40 iterations. 

To construct the halo we used again 50 iterations and 
the remaining parameters were taken as follows : ti = 50, 
dt = 1/2 4 , e = 0.04 and 9 t = 0.9. For fixing the veloc- 
ity anisotropy radial profile (|13[1 we used the algorithm de- 
scribed in section 12.3.21 The number of layers in this algo- 
rithm was ndiv = 500. 

Once all three components of our model were con- 
structed, we simply stacked them in order to obtain the 
complete system. In order to check whether this was indeed 
near equilibrium, as it should, we simply evolved with a full 
JV-body simulation, using again gyrfalcON, now with an in- 
tegration step and softening length of dt — 1/2 7 and e = 0.02 
(parameters were chos e n acc ording to recommendations of 
iRodionov fc Sotnikoval ((2005)). The tolerance parameter for 
gyrfalcON was set 8t = 0.6. The evolution of the total sys- 
tem over 160 time units is illustrated separately for the disk, 
bulge and halo components in figures[|)][7]and[S] respectively. 
These show that all three components of the constructed 
model conserve their structural and dynamical properties 
very well, demonstrating that the constructed model is in- 
deed close to equilibrium as it should. 

An interesting question arises in connection with the 
equilibrium of our model: how well do the moments of the 
velocity distribution in the c onstructed disk satisfy the equi- 
librium Jeans equations (see lBinnev fc Tremainelll987h ? 
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X 

Figure 10. The edge-on line-of-sight mean velocity v\ m (x) and 
and velocity dispersion O"i os (a;) for model DISK.SVR. These pa- 
rameters are defined in section 12.3.31 
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Here <l?tot is the potential generated by all the components 
of our model (disk, halo and bulge), v v , <tr, a v and a z are 
four moments of the velocity distribution in the disk (mean 
azimuthal velocity and velocity dispersions in the R, tp and 
z directions, respectively). 

Fig. [9] comparises the radial profiles of v v , <r v , and a z 
calculated from the constructed disk and from the Jeans 
equations (|14p . It can be seen that the model follows the 
Jeans equations very well. Note that the moments of veloc- 
ity distribution both in the Jeans equations and in the con- 
structed disk depend on z. In fig.[U all moments calculated 
by means of Jeans equations were calculated for z — 0. We 
thus took only particles with z < 0.05 to calculate them 
from simulations. We also checked that our disk follows the 
Jeans equations very well in the rest part of the space (not 
shown here). 

We want to underline that, according the Jeans equa- 
tions, the a z (R,z) in our disk is unambiguously defined by 
the chosen mass model of the galax y (see third equation of 
system (O), as already discussed bv lRodionov fc Sotnikoval 
(20061 1. 

3.3 Models with given line-of-sight kinematics 

In this section we demonstrate the capability of the iter- 
ative method to construct models with given line-of-sight 
kinematics. Let us first calculate the edge-on line-of-sight 
mean velocity v\ oa (x) and the edge-on line-of-sight velocity 
dispersion o"i os (x) of the disk model constructed in the previ- 
ous section (model DISK.SVR in section [3~2)l . These profiles 
are presented on figure [lOl 

We construct two disk models. The first one, called 
DISK.MVLOS, with a given v\ OB (x) and the second one, 
called DISK.SVLOS. with a given ai os (a;). Since we use the 
disk galaxy mass model described in the previous section, 
the process is similar to reconstructing DISK.SVR by using 
line-of-sight kinematic profiles obtained from "observation" . 
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Figure 6. Initial evolutionary stages for the disk of the constructed disk galaxy model. The evolution of the model was calculated with 
live disk, halo, and bulge components (see also fig. [7] and |8j. From left to right, the upper snapshots show the disc views face-on for 
times 0, 40, 80, 120 and 160 and the grey scales are logarithmically spaced. The middle and bottom panels show the dependence of 
various disc quantities on the cylindrical radius R at the same times. Her e n is the number of particle s in concentric cylindrical layers; 
z i/2 i s the median of the value \z\, i.e a measure of the disc thickness ( see ISotnikova &c Rodionovll2006T) and v^, o^j, and u z are four 
moments of the velocity distribution. At the beginning of the evolution (t = 0) the disk has, by construction, the radial dispersion profile 
given by eq. 11121 . 



Our initial model for the iterative method was a "cold" 
disk where all particles move on circular orbits (see previous 
section where we constructed DISK.SVR). We made 100 it- 
erations, each with ti — 20. Note that this is twice the num- 
ber of iterations used for DISK.SVR, because the converge 
in the case of line-of-sight kinematics is slower. The remain- 
ing parameters were taken dt = 1/2 4 , e = 0.04 and 6 t = 0.9. 
In order to fix the profile of v los (x) (for DISK.MVLOS) and 
the profile of <7i oa (:r;) (for DISK.SVLOS) we used the algo- 
rithms described in section 12.3.31 The number of layers in 
these algorithms was ndi v = 200. 

Let us check the equilibrium of the DISK.MVLOS and 
the DISK.SVLOS disks. In both cases we use the halo and 
bulge constructed in the previous section, because the mass 
model is the same. For the self-consistent evolution we used 
the same parameters as in previous section. The evolution of 
the disks in these models is illustrated in figures nTl and [T2l 
respectively. These show that the constructed disks are as 
close to equilibrium as they should. 



It is interesting to compare DISK.SVR and the two 
disks constructed from it by using line-of-sight kinematics. 
The three radial profiles of or and a v are visibly differ- 
ent, as can be seen in figure 1131 The velocity dispersion 
in the disk plane is visibly bigger for DISK.MVLOS than 
for DISK.SVR, especially near the disk periphery. This is 
also the case for DISK.SVLOS, but to a lesser extent. In 
general, it is clear that both models DISK.MVLOS and 
DISK.SVLOS are different from DISK.SVR. This was not 
expected and must be due to the fact that more than one 
equilibrium solution exists for the adopted constraints. This 
must be kept in mind when we apply our method for con- 
structing phase models of real galaxies and we will discuss 
it more extensively in a forthcoming paper. 
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Figure 7. Initial evolutionary stages for the bulge of the constructed disk galaxy model. The evolution of the model was calculated 
with live disk, halo, and bulge (see also fig. [6] and |8j. The upper snapshots show the bulge viewed edge-on for two moments of time 
(0, 160); the grey intensities correspond to the logarithms of particle numbers in the pixels. The bottom panels show the dependence 
of two parameters of the bulge on the spherical radius r for various moments of time. Here n is the number of particles in concentric 
spherical layers; a T is the velocity dispersion in the r direction (in spherical coordinate system). We demonstrate these parameters in 
order to show that the model is close to the equilibrium. For astrophysical applications it should be taken into account that the bulge 
in our model is not spherical. 



4 CONCLUSIONS 

We presented a new method for constructing equilibrium 
phase models for stellar systems — the iterative method. 
The aim of this method is to construct equilibrium TV-body 
models with given parameters, or constraints. More specif- 
ically, these are a given mass distribution and, if desired, 
given kinematic properties, parameters, or constraints. Our 
method is straightforward both conceptually and in its im- 
plementation. We believe that it is this simplicity that makes 
this method so powerful. It simply relies on a constrained, or 
guided evolution. We let the system reach equilibrium via 
a dynamical evolution in a number of successive steps. In 
between two such steps we make sure that the parameters 
are set to their desired value and/or that the constraints are 
fulfilled. This means that the evolution is guided towards an 
equilibrium with the desired parameters and/or constraints. 



Setting a mass distribution is of course obligatory, but kine- 
matical constraints are not. If we wish to include them, we 
have the choice of a large number of possibilities, such as 
setting the radial profile(s) of one, or more moments of the 
velocity distribution. In this article we described only a few 
types of kinematic constraints: the profile of radial velocity 
dispersion, the profile of velocity anisotropy, a condition of 
velocity isotropy and line-of-sight kinematics. Procedures for 
further types of kinematic constraints can be easily found 
following similar techniques. Furthermore, our implemen- 
tation of the iterative method can be directly applied to 
systems with arbitrary geometry, i.e. the given mass distri- 
bution can be arbitrary and need not have any symmetries. 
Thus our method can be used in many different applications. 

We used our iterative method to construct several mod- 
els. The first one is a triaxial system. The second one is 
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Figure 8. Initial evolutionary stages for the halo of the constructed disk galaxy model. The evolution of the model was calculated with 
live disk, halo, and bulge components (see also fig. [6] and 0. We show the dependence of various halo parameters on the spherical radius 
r for various moments of time. Here n is the number of particles in concentric spherical layers (upper left panel); crg/cr r is the ratio of 
velocity dispersion in the 9 and r directions (upper right panel); o> is velocity dispersion in the radial direction (bottom left panel); a v 
is the velocity dispersion in the tp direction (bottom right panel). At the initial moment of time the halo has the profile of ag/cr r given 
by 



a multi-component model of a disk galaxy consisting of a 
stellar disk with a given radial velocity dispersion profile, a 
non-spherical bulge and a halo with a given anisotropy pro- 
file. We also constructed two disk models with given line-of- 
sight kinematic. Using self-consistent A-body simulations, 
we made sure that the models we constructed are indeed 
very close to equilibrium (see figs. [3] [6] [7J [8] [12] and [XT) . 



The iterative method has a number of further applica- 
tions. It can of course be used for constructing equilibrium 
initial conditions for A-body modelling of stellar systems. 
For instance, the iterative method allows one to investigate 
bar formation in galaxies with a halo having different kine- 
matics. Also the iterative method can be applied for con- 
structing phase models of real galaxies. For example, we 
can model observational data by constructing phase models 
with given line-of-sight kinematics, as shown in section [3T3T 
This paves the way for studies of e.g. the distribution of dark 
matter in ellipticals, or obtaining phase space models of ob- 
served disk galaxies. A further interesting application is the 



study of the properties of several equilibrium models for a 
given mass distribution, as for example triaxial systems. 

The software necessary for the implementation of 
this method should be thought of in a very modular 
way, e.g. with different units for the various kinemati- 
cal constraints, and is very straightforward to write. Nev- 
ertheless, we will make our own software publicly avail- 
able as soon as this paper is accepted, at the address 



http://www.astro.spbu.ru/staff/seger/soft/. This package 



will contain also step-by-step examples for constructing 
models by using the iterative method, including the mod- 
els described in this article. Our sof tware uses the A-body 
code gyrfalcON l|Dehnenll200ol . 120021) and the N EMO pack- 
age | |http://astro.udm.edu/nemo| Teubenlll995l '). 
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Figure 9. Comparison of profiles of the velocity distribution moments calculated from the Jeans equations and from the disk of the 
constructed disk galaxy model (DISK.SVR). All moments for the disc were calculated inside the region \z\ < 0.05. Left panel: the solid 
line corresponds to the value v v for the model, and the dashed line corresponds to the same value calculated from the Jeans equation 
(the first equation of the system ||14|| ) for z = 0, and where the values an and a<p were taken from the model. Middle panel: the solid line 
corresponds to the value a v for the model, and the dashed line corresponds to the same value calculated from the Jeans equation (the 
second equation of the system JH}), where the values v v and or were taken from the model. Right panel: the solid line corresponds to 
the value a z for the model, and the dashed line corresponds to the same value calculated from the Jeans equation (the third equation of 
the system JH}). 
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Figure 13. The dependence of four moments of the velocity distribution, namely v v , otj, a v and c z as a function of the cylindrical 
radius R for DISK.SVR, DISK.MVLOS and DISK.SVLOS. 



